6.1.擬似相関_1

ゲーム時間とテスト成績の散布図

data <- read.csv("data/Ch6.data.csv")
plot(data[,c(2,1)], xlab="ゲーム時間", ylab="テスト成績")
lm(grade~hours, data) %>% abline(col=4)

相関係数は-0.49です

6.1.擬似相関_2

図6.2.を再現するには「家庭環境」データが必要なため、データ生成からやり直す

plot(data[,c(2,1)], xlab="ゲーム時間", ylab="テスト成績")
lm(grade~hours, data = data0) %>% abline(col=4)

相関係数は-0.45。 テキストとちょっと違うが、許容範囲かと。

6.1.擬似相関_3

3変数間の散布図。

library(psych) #散布図用
data0 %>% pairs.panels(rug=F,ellipses=F,lm=T)

6.1.P119_図6.2的な無向グラフ

library(qgraph) #無向グラフ用
data0 %>% cor() %>% qgraph(edge.labels=T)

6.1.3.同時性

pc <- read.csv("data/police_crime.csv", stringsAsFactors=FALSE)
plot(pc[,2],pc[,3],type="n", xlim = c(1.5, 3.5),
 xlab="警察官数", ylab="刑法犯認知件数")
text(pc[,2],pc[,3], pc[,1])
lm(crime~police, pc) %>% abline(col=4)

相関係数は0.129。

6.2.2.平均トリートメント効果_1

set.seed(1234); n <- 400; T <- rbinom(n, 1, 0.6)
TE <- 2; Y <- TE*T + rnorm(n); id <- 1:400
data1 <- data.frame(id, T, Y)
data1 %>% tidyr::spread(T, Y) %>% head() %>% kable()
id 0 1
1 NA 2.485227
2 0.6967688 NA
3 0.1855139 NA
4 0.7007335 NA
5 0.3116810 NA
6 0.7604624 NA

NAの部分の値がわかればトリートメント効果が計算できるが実際はわからない。

6.2.2.平均トリートメント効果_2

boxplot(Y ~ T, data=data1)

EY1 <- mean(Y[T==1]); EY0 <- mean(Y[T==0])

平均トリートメント効果(ATE)は1.872で実際の値(2)に近い。

6.2.2.平均トリートメント効果_3

set.seed(12345); n <- 400; Z <- runif(n)
T <- rbinom(n, 1, Z); TE <- 2; Y <- TE*T + (2*Z - 1) + rnorm(n)
data2 <- data.frame(T, Y, Z)
boxplot(Y ~ T, data=data2)

EY1 <- mean(Y[T==1]); EY0 <- mean(Y[T==0])

平均トリートメント効果(ATE)は2.813で実際の値(2)から外れている。

7.1_外生変数

data("ToothGrowth")
summary(ToothGrowth)
##       len        supp         dose      
##  Min.   : 4.20   OJ:30   Min.   :0.500  
##  1st Qu.:13.07   VC:30   1st Qu.:0.500  
##  Median :19.25           Median :1.000  
##  Mean   :18.81           Mean   :1.167  
##  3rd Qu.:25.27           3rd Qu.:2.000  
##  Max.   :33.90           Max.   :2.000

OJとVCそれぞれ30レコードずつ。 それぞれ、doseとlenの関係は?

ToothGrowth %>% lm(len~dose, data = .) -> lm.TG
ToothGrowth %>% filter(supp == "VC") %>% lm(len~dose, data = .) -> lm.VC
ToothGrowth %>% filter(supp == "OJ") %>% lm(len~dose, data = .) -> lm.OJ
summary(lm.TG)$coef
##             Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 7.422500  1.2600826  5.890487 2.064211e-07
## dose        9.763571  0.9525329 10.250114 1.232698e-14
summary(lm.VC)$coef
##             Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  3.29500   1.427060  2.308943 2.854201e-02
## dose        11.71571   1.078756 10.860392 1.509369e-11
summary(lm.OJ)$coef
##              Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 11.550000   1.721951 6.707508 2.788784e-07
## dose         7.811429   1.301673 6.001070 1.824801e-06
plot(ToothGrowth[,c(3, 1)], col = ToothGrowth$supp)
abline(lm.TG, col=4)
abline(lm.VC, col=2)
abline(lm.OJ)

VCはOJと比べて、doseが少ないときのlenの値が小さい。 doseが2を超えると、VCとOJの差は殆ど無い。

as.double(lm.TG$coef[1] + lm.TG$coef[2] *2.114)
## [1] 28.06269
as.double(lm.VC$coef[1] + lm.VC$coef[2] *2.114)
## [1] 28.06202
as.double(lm.OJ$coef[1] + lm.OJ$coef[2] *2.114)
## [1] 28.06336

7.2_内生変数

set.seed(1234)
n <- 200; e <- rnorm(n)
X <- (1+0.5*e)*runif(n)
b0 <- 1; b1 <- 2
Y <- b0 + X*b1 + e
lm(Y~X)$coef
## (Intercept)           X 
##  0.02967807  3.71966851

Xの係数3.7196685はb1の推定としては不適切?

cov(X,e)/var(X) + b1
## [1] 3.719669

数式(7.7)による補正rをすると、上記Xの係数とほぼ一致

7.3.2_双子データ

OpenMxパッケージの双子データを用いる

library(OpenMx)
data(twinData);head(twinData)
##   fam age zyg part wt1 wt2    ht1    ht2   htwt1   htwt2    bmi1    bmi2
## 1 115  21   1    2  58  57 1.7000 1.7000 20.0692 19.7232 20.9943 20.8726
## 2 121  24   1    2  54  53 1.6299 1.6299 20.3244 19.9481 21.0828 20.9519
## 3 158  21   1    2  55  50 1.6499 1.6799 20.2020 17.7154 21.0405 20.1210
## 4 172  21   1    2  66  76 1.5698 1.6499 26.7759 27.9155 23.0125 23.3043
## 5 182  19   1    2  50  48 1.6099 1.6299 19.2894 18.0662 20.7169 20.2583
## 6 199  26   1    2  60  60 1.5999 1.5698 23.4375 24.3418 22.0804 22.3454
##    cohort zygosity age1 age2
## 1 younger     MZFF   21   21
## 2 younger     MZFF   24   24
## 3 younger     MZFF   21   21
## 4 younger     MZFF   21   21
## 5 younger     MZFF   19   19
## 6 younger     MZFF   26   26

fam:家族番号
age:双子の年齢
zyg:遺伝子とコホートをコード化(1-10)
part:?(0:15, 1:224, 2:3569)
wt1,2:双子の体重(kg)
ht1,2:双子の身長(m)
htwt1,2:双子のBMI(kg/m^2)
bmi1,2:双子のBMIの対数
cohort:youngerかolder
zygosity:遺伝子パターン(MZFF, MZMM, DZFF, DZMM, DZOS)

boxplot(age~cohort, twinData)

cohortの境界は30歳

双子の組に対して、身長の差と体重の差で線形モデルを作る

twinData %>%
  mutate(wt0 = wt1 - wt2, ht0 = ht1 - ht2) -> data1
lm(wt0 ~ ht0-1, data=data1) -> lm.data1
summary(lm.data1)$coef
##     Estimate Std. Error  t value Pr(>|t|)
## ht0 87.56072   1.661656 52.69486        0
plot(data1$ht0, data1$wt0, xlab="身長", ylab="体重")
abline(lm.data1, col=4)

身長1cm毎に体重は0.87kg増える

c("age", "zyg", "part", "wt0", "ht0", "cohort", "zygosity") -> coln
rbind(
  twinData[,c(2:4,5,7,13:14)] %>% set_colnames(coln),
  twinData[,c(2:4,6,8,13:14)] %>% set_colnames(coln)
)-> data2
lm(wt0 ~ ., data=data2) -> lm.data2
summary(lm.data2)$coef
##                  Estimate Std. Error     t value      Pr(>|t|)
## (Intercept)   -67.0140567 2.42647491 -27.6178651 7.880592e-160
## age             0.1270048 0.01097231  11.5750260  1.017459e-30
## zyg             0.6412055 0.07149348   8.9687268  3.752456e-19
## part            0.1472285 0.57012365   0.2582397  7.962292e-01
## ht0            72.6191300 1.31619664  55.1734654  0.000000e+00
## cohortyounger   1.0833814 0.46770506   2.3163773  2.056505e-02
## zygosityMZMM    3.0815594 0.32938609   9.3554632  1.084290e-20
## zygosityDZFF   -1.0652820 0.26069244  -4.0863555  4.428643e-05
## zygosityDZMM    2.4107945 0.37732550   6.3891639  1.769711e-10

身長1cm毎に体重は0.72kg増える。